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Abstract. Exact and nonperturbative quantum master equation can be 
constructed via the calculus on path integral. It results in hierarchical equations 
of motion for the reduced density operator. Involved are also a set of well- 
defined auxiliary density operators that resolve not just system-bath coupling 
strength but also memory. In this work, we scale these auxiliary operators 
individually to achieve a uniform error tolerance, as set by the reduced density 
operator. An efficient propagator is then proposed to the hierarchical Liouville— 
space dynamics of quantum dissipation. Numerically exact studies are carried out 
on the dephasing effect on population transfer in the simple stimulated Raman 
adiabatic passage scheme. We also make assessments on several perturbative 
theories for their applicabilities in the present system of study. 



1. Introduction 

The central problem of quantum dissipation theory is to study the dynamics of 
quantum system embedded in quantum thermal bath. The primary quantity of 
interest here is the reduced density operator, p{t) = trBjOT(i), after the bath degrees 
of freedom are all traced out from the total composite density operator. Due to its 
fundamental importance, quantum dissipation theory has remained as an active topic 
in diversified fields [1-10]. The challenge here, from both formulation and numerical 
aspects, is nonperturbative dissipation, with multiple time scales of memory, under 
time-dependent external field driving. 

For Gaussian stochastic force, the infiuence of bath on system can be characterized 
by force-force correlation functions. Exact formalism had then been established via 
the Feynman- Vernon influence functional approach [1-5]. Direct numerical integration 
methods, based on discretization of the path integral and summation up of the memory 
correlated terms, have been put forward such as the quasi-adiabatic propagator 
method [11-14] or the real-time quantum Monte Carlo scheme [15-19]. 
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The alternative is the differential approach, especially in a linear form to maximize 
the numerical advantage. It has also the advantage in the study of various dynamics 
such as the spectroscopic or control problems [20]. The calculus-on-path-integral 
(COPI) method is hence proposed to construct the differential counterpart of the path 
integral theory, reported as the hierarchical equations of motion (HEOM) formalism 
[21-30]. This formalism can also be derived via the stochastic description of quantum 
dissipation [31-35]. The COPI algorithm provides a unified approach to the influence 
of quantum environment ensembles, either canonical or grand canonical, and either 
bosonic or fermionic [551 US]- The COPI algorithm also takes into account the 
combined effects of multiple memory time scales, system-bath coupling strengths, and 
system anharmonicity. The resulting HEOM formalism is therefore nonperturbative in 
nature, and always converges in principle. Moreover, the HEOM formalism is exact, 
not just for its propagation equivalent to the path integral theory, but also for the 
fact that the initial correlations between system and bath can now be incorporated by 
the steady-state solutions to the HEOM, before external time-dependent fields taking 
effect. Recently, we have further developed a numerical efficient filtering method for 
the propagation of the HEOM [361137]. 

In this work, we report a HEOM-based study on population transfer with 
dephasing in the scheme of stimulated Raman adiabatic passage (STIRAP) 38] . The 
laser control of dissipative systems has been addressed extensively [39-49] , but mostly 
on the basis of weak dissipation treatment. The correlated influence of driving and 
dissipation is often important, as demonstrated previously [50l|5T]- With the aid of 
the numerically exact results, we analyze the dephasing effects on transfer dynamics 
in relation to the STIRAP mechanism and examine some second-order quantum 
dissipation theories for their applicabilities in the systems of study. 

The remainder of paper is organized as follows. We present the HEOM formalism 
together with comments on its numerical implementation in Sec.[2l and the derivations 
in Appendix. In Sec.[3l we study the dephasing effect on population transfer dynamics 
in the STIRAP scheme. We report the numerically exact results via the HEOM 
formalism, followed by discussions in relation to the STIRAP mechanism. In Sec. HI 
we present the details of numerical performance of the HEOM results, and make 
concrete assessments on several approximated quantum dissipation theories. Finally 
we conclude the paper. 

2. Hierarchical equations of motion formalism for quantum dissipation 

2.1. Description of stochastic bath coupling 

The total system-plus-bath Hamiltonian can be written in general as 



The last term denotes the multi-mode system-bath interactions. The involving 
system operators {Qa} are called the dissipative modes, through which the generalized 
Langevin forces {Fa{t) — e*''^*i^ae~'''^*} from the bath {hg) act on the system. For 
convenience, let the dissipative modes be dimensionless. The time dependence in the 
system H{t) arises from external driving fields. Throughout this paper, we denote the 
inverse temperature (3 = l/^k-sT) and set h=l. 




(1) 



a 
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We treat the Langevin forces as Gaussian stochastic processes. Therefore their 
effects on the system are completely characterized by the correlation functions, 

Cabit-T) = {Fjt)F,{T))^. (2) 

Here, (0)b = i'^eiOp^) denotes the thermodynamics average over the canonical 
ensembles of the bosonic bath. The correlation functions satisfy the symmetry and 
detailed-balance relations, or equivalently the fluctuation-dissipation theorem [21110]: 

with Jab{^) = —Jbai—^) = ^bai^) being the bath spectral density functions. The 
HEOM formalism requires Cab{t) be expanded in certain series form, so that the 
hierarchy can be constructed via consecutive time derivatives on path integral. Various 
schemes [22l [281 Ell [53; have been proposed to expand Cab{t) in exponential series, 
on the basis of analytical continuation evaluation of Eq. In particular, the hybrid 
scheme that also exploits quadrature integration method is applicable for arbitrary 
spectral density functions (30j . 

For simplicity we set Cab{t) — Caa{t)Sab- In this case the contributions from 
different dissipative modes {Qa} are additive. Without loss of generality, we present 
the formalism explicitly only for the single-dissipative-mode case, Qa — Q. We thus 
omit the index a for clarity of formulation. We also adopt the super-Drude model, 

J{uj) = . ■„ . (4) 

[(^/7)2 + 1]2 



The corresponding correlation function can be analytically evaluated as [20l [28l [53] 

M 

C{t >0) = [v+ {Pr + 7i]e-^* + Y J>me-^'"* + SCit). (5) 



m— 1 



All coefficients here are real and given in Appendix [cf. Eqs. ((AJll and ([AlOl) ]. The first 
term arises from pole of the spectral density function, which is of rank two. The second 
term is from the Matsubara poles, with 7,„ = 27rm//3 being the Matsubara frequency. 
The last term is the Matsubara residue, which would approach to zero if M — > oo. In 
this work, we adopt the Markovian residue ansatz j25n3Sj, i.e., 7me~^'"* |„j>jv/ ^ ^^^)'^ 
thus, 

oo _ _ M _ 

5C{t)^A5{t); - = l-'^-Y.-- (6) 

, 1 7m P 7 1 7m 

m=M+l ' ' m=l ' 



2.2. The HEOM formalism 

The dynamics quantities in the HEOM formalism are the reduced density operator p{t) 
and a set of auxiliary density operators (ADOs), {p„{t)}, that hierarchically resolve the 
memory contents of the bath correlation functions in the exponential series expansion 
of Eq. ([5]) . The index n that specifies an iV*''-tier ADO p„ consists of a series of 
nonnegative integers, 

n = {n, n', n, n', ni, • • • , nA/} , with n + + n + n' + rii + • • • + fiM = N. (7) 
Comparing to the reduced density operator p(t) = po{t) of primary interest, the 
specified p„ would have the order of • \i'r + ' nm=i I'^ml"", for its 

dependence on the individual components of interaction bath correlation functions in 
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the series expansion of Eq. ([5]) . These scahng factors will be incorporated properly in 
the final dimensionless p^, in order to validate a filtering algorithm for the numerical 
efficiency of the HEOM formalism. On the other hand, the indices in the set n of 
Eq. ((T]) cover all accessible derivatives of the Feynman- Vernon influence functional; 
see Appendix for the details. 

The final HEOM formalism is summarized as follows. It has the generic form of 

= - [iC{t) + r, + 6TZ]p„ + pV' + pir' + . (8) 

Here, C{t)pn = [H{t),pn], which depends in general on the external driving fields; T„ 
is the damping parameter that collects all related exponents, and 6TZ is the residue 
dissipation superoperator due to SC (t) . For the bath correlation function in the series 
expansion, Eq. (O with Eq. ([5]), they are given respectively by 



Fn = (n + n' 



M 

E 



Snp, = A[Q,[Q,p,]]. 



(9) 



Apparently, Fq = Fn|Ar=o = 0. 

The last three terms in Eq. ([8]) denote how the specified iV*''-tier ADO pr, depends 
on other ADOs of the same tier, the (A^— l)*''-tier, and the (A^ +l)*^-tier, respectively. 
For the bath correlation function in Eq. ([5]) , they are given explicitly by 



PrT — KiPn + K'Pn' : 



M 

'i[Q,\nP„-] + {Q,K'Pn'-} - 'i[Q, E -^""^n;; 



(10) 



Q,Ki+iPu+ +sign(i/r)-A^+iPn+ - Xf,,+iPi^'+ 



M 



Xn„, + lPn+ 



Here, A„ 



\/ri\^\ , Afi,^ = ^/fh^^p^\, = y/n\U^\, and A^* = 7\i+i\/«7H'> 
with the italic-font indices being from those in n of Eq. ^ . The indexes variations in 
Eq. (|10[) that specify those ADOs participating in the equation of pn are exemplified 
as follows: 



n = {n + 1, n', n — 1, n', ni, • • • , tlm} , {n ± 1, n', n, n', hi, • • • , hm} ■ (H) 

Similarly, n' differs from n of Eq. ([7]) only by changing {n',n') to (ti' + 1, n' — 1), while 
by changing to fim ± 1, and so on. Also note that p^ is an Af*'^-tier ADO, 
while Pn± is of an {N ± 1)*'' tier, as inferred from the second identity of Eq. ([7]). 

The initial conditions to the HEOM in the study of driven dissipative dynamics 
are obtained via the steady-state solutions to Eq. ([5]), before the time-dependent 
external fields interactions. For the steady-state solutions satisfying = 0, Eq. ^ 
reduces to a set of linear equations, under the constraint of Trpo = 1- The resulting 
Pn* is used as the initial Pn(io) to the HEOM. The initial system-bath correlations are 
accounted for by those nonzero initial ADOs. 



2.3. Comments on numerical implementation 

For the numerical HEOM propagation, we would like to have certain convenient 
working index scheme to track the multiple indices, denoted now as an ordered set 
of n = {ni, • • • , uk}, that specifies pn- Here we will provide two such schemes. The 
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number of the A^^'^-tier ADOs, with ni + ■ ■ ■ + uk — N, is ^IJt^^_ly = [k]- 
scheme, the ADOs are arranged as pn = pj^ with jn initiahzed by jn=o = and then 

N-l K N-Sk 

jn^O = jm-UK ^ ^ [k']+^ ^ [K-k]'^ Sk = l + ni + --- + nk. (12) 

N'=0 fc=l q=0 

Let L be the maximum level of the hierarchical tier. The total number of the ADOs 
{p.;0<N <L} is 

-^-E[^-]-{i-}- (13) 
In another scheme, ADOs can also be arranged as pn = p;„; = 0, • • • ,A/' — 1, with 

= z„^...„^ = -1 + E E { -^"'^^ } • (14) 

k=2 q=l 

Both schemes, Eq. (IT2]) and Eq. p^ . allow easy tracking of the coupled ADOs in the 
HEOM. The former [Eq. (fT2)l ] is somewhat more convenient in the filtering propagator 
described soon since it does not depend on L. 

The major difficulty in implementing the HEOM formalism is its numerical 
tractability. The number of ADOs, Af of Eq. itself alone can be huge in the 
case of strong non-Mar kovian system-bath coupling and/or low temperature, as large 
L and/or large K implied. Thus, a brute-force implementation is greatly limited by 
the memory and central processing unit (CPU) capability of computer facility, even 
for a two-level system where each ADO is a 2 x 2 matrix. 

To facilitate this problem, Shi, Xu, Yan and coworkers have recently proposed 
an efficient numerical filtering algorithm that often reduces the effective number of 
ADOs by order of magnitude [36l [37] . In reality, there is usually only a very small 
fraction of total ADOs significant to the reduced system dynamics. To validate the 
accuracy-controlled numerical filtering algorithm, the present HEOM formalism has 
been scaled properly so that all ADOs {pnit)} are of a uniform error tolerance. This 
remarkable feature is suggested by comparing the HEOM theory with the stochastic 
bath interaction field approach in the case of Gaussian-Markovian dissipation |37j . 
The involving ADOs are just the expansion coefficients, over the normalized harmonic 
wave functions that are used as the basis set for resolving the diffusive bath field [37] . 
Our numerical HEOM propagator exploits the filtering algorithm [36' . It goes simply 
as follows. If a pn{t) whose matrix elements amplitudes become all smaller than the 
pre-chosen error tolerance, it is set to be zero. Apparently, the filtering algorithm 
also automatically truncate the required hierarchy level on-the-fly during numerical 
propagation. By far the truncation for the Matsubara expansion still goes by checking 
convergency. 

3. Effect of dephasing on population transfer via STIRAP 

3.1. Numerical results 

The STIRAP is celebrated as an efficient and robust method for population transfer 
[5H] . It is characterized by its counterintuitive field configuration. For a three- level 
A-system as Fig.[Tl the Stokes pulse proceeds the pump pulse, and the intermediate 
state remains effective in dark. The STIRAP mechanism [38] is rooted at the existence 
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Figure 1. (a) A schematic view of the STIRAP of a three-level A system, (b) 
Population transfer under the STIRAP scheme for the dissipation-free gas phase. 
The time is in the unit of /3. See parameters in the text. 



of the coherent population trapping state under the two-photon resonance condition, 
ivg — LUp = ei — 63, in the A-system. Dephasing destroys this condition in terms 
of resonance and/or the existence of coherent population trapping state. The effect 
of dephasing on the simple STIRAP scheme has been studied extensively, but with 
approximations. These include phenomenological/perturbative methods [44-46], or 
classical/stochastical bath treatments [47-49]. 

We revisit the dephasing effect on the simple STIRAP-based population transfer, 
as the exact dissipative dynamics are now established with the present HEOM 
formalism. We also examine three schemes of second-order approximation [20, 52-55]: 

(i) The Redfield theory, which neglects the correlated driving-and-dissipation effect; 

(ii) CS-COP, which is the conventional time-nonlocal quantum master equation, 
including the field-dressed dissipation contribution, and equivalent to the present 
HEOM truncated at the first tier; (iii) CODDE, in which the driving field-free part of 
dissipation superoperator is time-local, while the field-dressed part is time-nonlocal. 
Neglecting the latter leads it to the Redfield theory. 

The total Hamiltonian under the rotating wave approximation assumes 

Hr{t) npit)Dp + ns{t)Ds + E + ~ E ^) 1 • (15) 



Here, Dp = |1)(2| + |2)(1| and Ds = |2)(3| + |3)(2|, while 17p(t) and flsit) denote the 
Rabi frequencies of the resonant pump and Stokes fields, respectively. The dissipative 
mode Qa = is responsible for dephasing. The interaction spectral density 

function Ja(w) = (7r/2) [c^j/(mjWj)]^(a; — ujj) assumes super-Drude as Eq. ([4]). 
The system Hamiltonian is then 

3 

H{t) = np{t)Dp + ns{t)Ds +Y,6ea\a){a\ . (16) 

The Caldeira-Leggett renormalization energy [56l [57] is Sea — ^ f^^du Ja{uj) / uj = 
jVala, for the super-Drude model [Eq. Q]. In the STIRAP configuration, it would 
relate to the effective detuning at short-time of the pump or Stokes field, as inferred 
from the analytical result of driven Brownian oscillator [5D1 [55] . 
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Figure 2. Evolutions of P33 and p22 (in insets) via the exact HEOM (solid), 
CODDE (dash), CS-COP (dot) and the Redfield equation (thin-solid) for single- 
dissipative— mode case: (a) Qi = |1)(1| and (b) Q2 = |2)(2|. The system— bath 
coupling strength r] = 0.64. The parameter /B-y = 5 exemplifies the Markovian 
condition. The time is in the unit of /3. 
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Figure 3. Same as Fig. [2] except the parameter f^j = 0.5 exemplifies the non- 
Markovian condition. The time is in the unit of f3. 



We set the pump and Stokes fields be of same Gaussian shape, ilp(t + tp) — 
0,s{t + ts) = Aexp[-^{wtf], but center them at tp = 200/3 and ts = -200/3, 
respectively and counter-intuitively. The driving strength and inverse duration 
parameters are set to be (3A = 0.1 and f3w = 0.005. The corresponding dissipation- 
free transfer dynamics is shown in Fig.[ljb). As here the bath influence is considered 
to be pure-dephasing in the absence of fields, the initial system is just chosen to be 
completely on the |1) state and all the ADOs are zero. For the effect of bath, we 
set the coupling strength rj — 0.64 [cf. Eq. (jl])], and consider both the Markovian and 
non-Markovian cases as follows. 

The Markovian transfer dynamics, under the influence of single dephasing mode 
of either Qi — or Q2 = |2)(2|, is exemplified in Fig.dl with /37 — 5. We observe: 

(i) The Qi-modc effect shown in Fig.[2Ua) leads to all three populations about 1/3 
after the driving; (ii) The Q2-inode effect shown in Fig.[3fb) is less sensitive than its 
Q\ counterpart, achieving a higher transfer efficiency, despite it is only about 0.55. 

The non-Markovian transfer dynamics is exemplified in Fig.[3l with = 0.5. 
In comparison with the Markovian counterparts, we observe: (iii) The Qi-mode case 
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Figure 4. Two-mode dissipation cases; (a) Qi = |1)(1| and Q2 = |2)(2| together, 
and (b) Qi = |1)(1| and Q3 = |3)(3| together, evaluated with /37 = 0.5. The 
system-bath coupling strength -q = 0.64. The time is in the unit of /3. 



behaves about the same; (iv) But the Q2-inode results in a higher transfer yield, 
increasing to about 0.73 via the exact calculation. We have also calculated the 
influences of Q3 = |3)(3| for both Markovian and non-Markovian cases. The results 
(not shown here) are similar to those of Qi, except for some small oscillations. 

Two double-modes {Q1+Q2 and Q1+Q3 uncorrelated) non-Markovian dephasing 
dynamics are shown in Fig. Ufa) and (b), respectively. They are insensitive to the 
non-Markovian parameter, and both reach at the final equal-populations, based on 
the numerically exact results. Comments on the approximated schemes, the CODDE, 
CS~COP, and Redfield theory, presented in Figs. [2]-[4] will be given later; see Sec. 14. 21 

3.2. Discussions 

The above observations can be understood by the well-established STIRAP 
mechanism [38]. The Qi-mode, which associates with the fluctuation of level easily 
destroy the two-photon resonance (TPR) condition, as described at the beginning of 
Sec.[3l Thus, it ends up with the observed equal populations in all accessible levels 
by the strong fields, as consistent with the analysis in [33] ■ The similarity between 
the Qs and Qi influences is also explained. The same reason accounts further for the 
case of uncorrelated two modes {Qi + Q2 or Qi + Q3) dephasing, as depicted in Fig. [4] 
It is anticipated that when 7 <C w (termed as the linear adiabatic limit below), the 
equal-population will be broken to be in favor of |3), due to the marginally partial 
fulfilment of the TPR condition. 

On the other hand, the (52-mode is associated with the fluctuation of the 
intermediate level |2). It alone does not affect the TPR condition. However, this 
condition, based on the numerically exact results shown in this work, is not sufficient to 
retain the coherent population trapping state, chosen ad hoc earlier for the dephasing- 
free STIRAP scenario in Fig.[T][b). It is anticipated that the coherent population 
trapping state may be recovered in the aforementioned linear adiabatic limit. This 
is in line with the observation-(iv), where the non-Markovian population transfer 
with single (52-inode dephasing [Fig.[3ljb)] is of higher efficiency than its Markovian 
counterpart [Fig.[2l[b)] . The previous study based on perturbative dephasing dynamics 
[44] has also shown that the single (52-niode does not affect the transfer efficiency in the 
linear adiabatic limit. Nevertheless, STIRAP in the presence of complex dephasing. 
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Table 1. Performance of HEOM formalism with filtering. 
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CPU (min) 


ATmax (AA) 
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Fig. 2(b) 


932 


19765 (8.44 X 106) 


17 


0.95 


1.28 


Fig. 3(b) 


15 


400 (9.24 X 10") 


9 


1.24 


348 


Fig. 4(b) 


266 


3664 (1.00 X 10'') 


9 


1.24 


348 



if 100% transfer ever achievable, would require dynamics feedback control of pump 
or Stokes laser frequency This would involve chirp and realize STIRAP in a 

nonlinear adiabatic condition, rather than the linear simplification considered here. 

4. Assessments on theoretical methods and concluding remarks 
4.I. Numerical performance of the HEOM formalism 

The numerical performance of the HEOM formalism with filtering is summarized in 
Table [1] for the systems reported in the three figures' (b)-panels. The CPU time is 
for a single Intel(R) Xeon(R) processor@3.00GHz to calculate the exact result in each 
(b)-panel for the time period -1200/3 < t < 2000/3 with the time step dt = 0.01/3 
using the fourth-Order Runge-Kutta propagator; A/'max denotes the largest number of 
active ADOs and L^ax the highest tier level, ever survived in the entire time span of 
the numerical propagation. The filtering error tolerance is chosen to be 10^^, following 
our previous work [36] . We input M = 6 for the number of Matsubara terms being 
explicitly included, which has been tested to give converged results of p{t) ~ Po{t) in 
all calculations. The total number M of mathematical ADOs follows Eq. (|13|) and is 
given inside the parentheses. The effect of filtering is clearly seen. The number of 
active ADOs with filtering is insensitive to the input M, as long as it is large enough. 
In the present study, the number of active ADOs reaches A/'max only during the period 
about —250 (3 < t < 500 /3 and grows up or drops down dramatically outside that 
period with the fields turning on or getting over. Apparently, A/'max increases with the 
number of dissipative modes. 

At least one (Lmax)*^-tier ADO actively participates during the HEOM 
propagation. Its leading contribution to the reduced density operator is of (2Lmax)*'^ 
order in the system-bath interaction. Physically, imax is closely related to the 
modulation K-parameter [27], introduced originally by Kubo for motional narrowing 
problem [6]. This dimensionless parameter is determined via k = or similar, 

for each individual exponential component in Eq. ([5]). The last two columns of Table 
[1] are the modulation parameters k and Km=i of the leading Matsubara term. The 
modulation K-parameter relation to the value of Lmax |27| can be clearly seen. In 
both the Markovian and non-Mar kovian cases of the present study, km = 7m/ \J\vm\ 
monotonically increases with m, cf. Eq. (jA.10|l . Actually, the Matsubara series 
truncation M in Eq. ([S]) can be estimated via its reaching the fast modulation 
condition, ku 3> 1. As the temperature decreases, km getting smaller, and eventually 
cause the value of Lmax be pretty large. The present HEOM construction is based on 
the Matsubara series expansion, which may no longer be numerically implementable 
in the extremely low temperature regime. Alternative expansion method such as the 
hybrid scheme [3D] is needed to the required HEOM construction. 
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4-2. Assessments on three second-order approximated theories 

With the exact results, we can now make concrete assessments on the three second- 
order approximated schemes, the Redfield theory, CS-COP, and CODDE, exploited 
in the numerical demonstrations. The dissipative modes {Qa — \a){a\} considered in 
Sec. [3] are all of pure dephasing, in the absence of external fields. The Redfield theory 
would be exact if there were no correlated driving-and-dissipation effect [201 l54] . 
Therefore, the non-Mar kovian dynamics manifest here as the correlated driving and 
dissipation. Apparently, the Redfield theory, by its Markovian nature, is independent 
of the width 7-parameter of bath spectral density. Observed is also the fact that the 
schemes of approximation are sensitive to the Q2-niode rather than the Qi- or Q3- 
mode dephasing. This fact is also easily understood, by considering their relations to 
the STIRAP mechanism as discussed earlier. Further remarks on the approximated 
theories for their applicabilities in the systems of study are as follows. 

The CS-COP theory [20, 531[5S] is overall most unsatisfactory, despite it contains 
formally a description of memory and driving-and-dissipation correlation. Even in 
the non-Markovian (52-niode case of Fig.[3l[b) where it appears to be superior than 
the Redfield theory, the CS-COP results in a decreased transfer efficiency from its 
Markovian counterpart shown in Fig.[2jb). This is qualitatively contradictory to the 
physical anticipation, as discussed earlier. 

The CODDE [20,,53ji55j appears to be the most favorable perturbation theory. It 
gives the best approximated transfer dynamics in all cases presented in Sec.[3l except 
the one to be discussed soon. Its overall superiority is also true in the driven Brownian 
oscillator systems [20l [51]. The CODDE is actually a modified Redfield theory, 
with inclusion of correlated driving-and-dissipation effects. The involving field- 
dressed dissipation kernel is time-nonlocal but constructed with a partial ordering 
resummation, rather than the chronological ordering prescription that characterizes 
the CS-COP [Soma [55]. 

The only exception is the Markovian Qi-mode case shown in Fig.[2l^a), where the 
Redfield dynamics is almost exact. The reason for this exception is also accountable. 
As we mentioned earlier, the non-Markovian dynamics manifest as the correlated 
driving and dissipation. This correlated effect diminishes in both the fast- and 
slow-modulation regimes, as inferred from the exact and analytical results of driven 
Brownian oscillator systems [58]. This conclusion can be carried over to the present 
system of study, as suggested here. Apparently, the identical value of = 5, adopted 
in the two cases of Fig. [21 acquires the fast-modulation limit for the Qi-mode, but 
not yet for the Q2-niode. In the latter case, the CODDE resumes its superiority. 

4.3. Closing remarks 

In summary, we have presented a hierarchical Liouville-space approach, which is 
exact and also quite tractable numerically, to general quantum dissipation systems 
under external driving fields. The auxiliary density operators are all of a uniform 
error tolerance, as that of the reduced density operator. We also make comments 
on the numerical facilitation about the multiple-index assignment and the filtering 
algorithm. We numerically study the dephasing effects on the population transfer, 
with a fixed simple STIRAP configuration, and present a concrete assessment on 
various approximation schemes. 
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Appendix. Construction of HEOM via the COPI approach 

This appendix gives the details of the COPI approach to the HEOM formalism. It 
starts with the influence functional in path integral. Let {\a)} be a basis set in the 
system subspace and set a = {a, a') for abbreviation. Denote the evolution of reduced 
system density operator in the a-representation by 

p{a,t) = p{a,a',t) = J daQU{a,t; ao,to)piao,to)- (A.l) 
Here, the reduced Liouville-space propagator is 

U{cx,t;aQ,to)^ / Pa e*^["lj^[Q;]e-*'^[" (A.2) 
Q;o[to] 

The effects of bath on the reduced system are contained completely in the influence 
functional T. For Gaussian stochastic forces {Fa{t)} from fluctuating bath, it assumes 
the Feynman- Vernon form [T], which can be recast as [27-29]: 



T[a]=cxp{- I dT^A[a(T)]Sa(T;{a})}, (A.3) 



with 



Aa[oc{t)] = Qa[a{t)] - QaWit)], (A.4) 
Bait; {a}) = Bait; {a}) - B'^it; {a'}), (A.5) 



and 



Bait;{a]) = Y. I dTCabit-T)Qk[aiT)l 
b 

B'ait;{a'})^Yl I dTC:bit^T)Qb[a'iT)]. 



b 



(A.6) 



to 



Here, Cabit) is the bath correlation functions, defined by Eq. ([2]). The functional 
Aa [Eq. (|A.4|) ] depends only on the local time and its operator level form is just the 
commutator of the dissipative mode Qa- 

The functional Ba [Eq. (|A.5p ] does however contain memory, which is to be 
resolved via the COPI algorithm of consecutive time derivatives on all memory- 
contained functionals. To construct a close set of HEOM via the COPI algebra, 
it needs a proper expansion of Cabit) such as the exponential series, while maintaining 
the fluctuation-dissipation theorem of Eq. ([3|) . A super-Drude parametrization scheme 
and the resulted HEOM formalism have been presented in our previous work |28j . 
A hybrid scheme, exploiting the analytical continuation and quadrature integration 
methods to evaluate the fluctuation-dissipation theorem, has also been proposed [3D] . 

To illustrate the COPI algorithm, consider the super-Drude model, 

l]abUJ + iT(ab<J^ 

Jabi^) - TT—, ^2 I 112 ■ \^-^> 
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The parameters are all real and satisfy the symmetry relations of rjat — "Hba^ lab = 76a j 
and ?7^j — — as inferred from the symmetric relations implied in Jati^^)- The 
resulting correlation functions via Eq. ([3]) are 

M 

Cab{t > 0) = [{i^f+tiyf) + {Df+tDf)^abt]e-^^'' + Y, KnC-^^-' +6Cab{t).{A.8) 

m— 1 

The second term arises from the Matsubara poles, with 7m = 27rm//3 the Matsubara 
frequencies, and — —i{2/ (3)Jab{—ilm) = {^m)* is real, as inferred from the 
symmetry relation of the spectral density function in analytical continuation. The 
first term arises from pole of the spectral density function of rank 2. The involving 
coefficients are summarized as follows [20l HH [53] . 

= V'abllbl^. - -\{riab + V'ablabhlb , = "^f C0i{(3^ ab / 2) , 

yf- = -yf COt(/37,b/2) - i>f (/37a6/2) CSc2(/37,fe/2) , 

and 

^ 2(7?afc7m + Vla^ll,) _ .2 -gfc gb ^ fab + ^Ib7m)/(m7r) 

/3[(7™/7a6)2 - 1]2 - ' ^™ [(7m/7ab)2 - 1]^ ' ^ ' 

The residue SCab{t) can be approximated via the Markovian ansatz, i.e., dCab{t) — 
AabS{t); cf. Eq.®. 

To proceed we denote for every distinct exponent terms in Eq. ()A.8p . 



Bab{t;{a})= / dTe-^-(*-^)Qb[a(T)], 

Jto 

Babit;{a})= [ dT-/ab{t-T)e-^-'^'-^'^Qb[a{T)l (A.ll) 

Jto 

B:^{t;{a}) ^ J^C' rdTe-^'"(*--)gb[a(r)]; m = l,---,M. 
b -^to 

They are related to the influence generating functionals as [cf. Eq. (|A.6p ] 
Bab= ~i{Bab~ B'^i,), B'ai, = Bab + B'^i,, 

Bab^-l{Bab~Kb). B',, = Bab + B'ab, (A.12) 

B:^^-z{B^~B':,); m = l,...,A/. 

Now we define the auxiliary density operators (ADOs) via 

p,{t)=Un{t,to)p{to), (A.13) 

where 

rait] 

'r,{a,t;ao,to)= / Vae''^^°'^T„[a]e-''^^°'\ (A.14) 

J an\tn] 



with 



= s^[\l[{BabT-^{B'aX'^^{Babf-^{B'a,)<^] \{{Bl,f-]T. (A.15) 



The scaling factor Sn is defined as 
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The ADOs {pn} defined in this way are dimensionless and possess uniform error 
tolerance to support the filtering method, as described in Sec. 12. 31 Since i]'^^ = 
hence = 0, we set = for the scaling factor Sn of Eq. (jA.16[) and hereafter. 
The index set for this complex multi-dissipative modes case is now 



n = {uab, n'^b^ Uab, n'^b^ fi^,---, J . 



(A.17) 



The HEOM can now be derived by taking the time derivative on p^,. Note that 
the COPI method is not just a time derivative technique, but provide a way to resolve 
the bath memory effect of the influence functional in the operative level. The final 
results are 



{-} 



where 



a.b a,m 

snd = Y,^ab[Qa,[Qb,o]]. 



(A.18) 

(A.19) 
(A.20) 



The details of the swap term plT^ , the tier-down p^ ^ and tier-up pj,^^ terms are 



{-} 



Pn 



= E V ("-f + ^)nab\i^?''/K''\ Pn<.. + V «6 + iXJ'^fVi^f I Pn'^, , (A.21) 



a.b 



Pn 



J2 \f^^\ [Qb , Pn:J + E \/<bWf\ {Qb , Pn- } 



M 



■ E E VK^ 7- E ^™ ' 



(A.22) 



Pn 



-*El'^"V(""^ + ^)/l''"'l Q''^P<. +Pf,J{nab + l)/\i^f\ Qa,P-„+^ 
a,b 



+i^fyiKb + ^)/\^i\ Qa,P-n'^ 



^T.<''y«b + ^)/\K''\ Qa,Pn 



a^b 



M 



(A.23) 



-*EE V^'C + 1 7m Q 

a m— 1 

The index variations here are similar to those described in Eq. (jlip . which is just the 
single-mode simplification. Equation (jlOp is thus obtained readily. 
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